In this chapter we will describe a method of stabilizing the advection-diffusion equations called Streamline Upwind Petrov-Galerkin (SUPG), introduced by prof. T. J. R. Hughes. This method is described, for example, in the work [1].
Let us recall first the advection-diffusion problem, from the example of the Eriksson-Johnson model problem, described, for example in [2].
This problem is defined in a square area \( \Omega = [0,1]^2 \) in the following way: we seek a function \( u \) such that
\( a(u,v)=l(v) \forall v \) where
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\int_{\Omega} \epsilon \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\int_{\Omega} \epsilon \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) is an extension of the Dirichlet boundary condition to the entire area, while \( \beta = (1,0) \) represents the wind blowing from left to right, while \( \epsilon = 10^{-6} \) is the diffusion coefficient.
The SUPG stabilization method (Streamline Upwind Petrov-Galerkin) modifies our equation by adding certain terms that do not change the "sense" of the equation, but enforce the inf-sup condition \( a(u,v) +{ \color{red}{\sum_K (R(u)+f,\tau \beta\cdot \nabla v)_K } }=l(v)+ {\color{blue}{\sum_K (f,\tau \beta\cdot \nabla v)_K } } \quad \forall v\in V \)
where
\( R(u)=\beta \cdot \nabla u + \epsilon \Delta u-f =\frac{\partial u}{\partial x }+\epsilon \Delta u-f \) represents the residual of our solution, that is, the error between the left and right-hand sides resulting from the discretization of our problem. In other words, on a continuous level, in the abstract infinite dimensional functional spaces (which theoretically allow for an exact approximation of the solution) the residual is zero (in abstract infinite dimensional spaces we do not go wrong). However, since we are using a computational grid and a finite number of approximation functions, a discrete-level residuum will measure us the error of the method.
\( \tau^{-1}=\left(\frac{\beta_x}{h_x} + \frac{\beta_y}{h_y} \right) + 3\epsilon \frac{1}{h_x^2+h_y^2} \)
\( \epsilon=10^{-6} \)
\( \beta = (1,0) \)
\( h_x,h_y \) is the element dimension (we divide our integrals into individual elements and on each element we calculate the modified integral using the formula containing the element diameter), and
\( (u,v)_K = \int_K uvdx \) denotes the dot product over the element \( K \).